Research on rheological behavior of fresh concrete single-cylinder pumping based on SPH-DEM

In contrast to traditional approaches to simulating fresh concrete, the model applied here allows issues such as liquid phase and the motion of sub-scale particles to be considered. The rheological behavior of fresh concrete materials was investigated, and the slump test and pumping process of fresh concrete were simulated by combining the smooth particle hydrodynamics coupled with discrete element method. Based on Bi-viscosity model and Bingham model, linear and nonlinear fitting of rheometer data and the derivation equations were educing. Bi-viscosity model and the Bingham model were compared in slump test. The results show that the Bi-viscosity model is more accurate in simulation, and the error percentage is less than 10%. The Bi-viscosity model was used to simulate and predict the results of slump experiment, and the influence of rheological parameters on the slump velocity and shape was obtained. The simulation analysis model of concrete single-cylinder pumping is established, and the experimental and simulation analysis models are compared. The results show that the SPH-DEM pumping pressure prediction is very close to the experimental results.

the pressure of the concrete pumping process, the results are one pressure value higher than the experimental values, which is explained by the fact that concrete undergoes compositional changes during transport in the pumping pipeline, which results in much smaller rheological parameters near the walls than the actual concrete.Secrieru et al. 18 and others, using the commercial software Fluent (Ansys Ltd), divided the concrete into two phases, using the rheological parameters of the lubrication layer near the wall and the internal plunger flow using the rheological parameters of the concrete itself for pumping pressure analysis, and the simulation results were very close to the experimental numerical results.
The phenomenon of formation and analysis out of the pumping process can also be used Discrete element method (DEM), which discrete the concrete into two phases of aggregate, slurry, and fitted with different material friction and recovery coefficients, respectively.Cui et al. 19 , discrete concrete large scale particles and continuous fluid, analyzed the concrete coarse aggregate shape for concrete test L box clogging factors were explored.Krenzer et al. 20 divided concrete into wet and dry particles, and applied particle replacement method to form a new material particle after wet and dry collision to simulate the aggregation of slurry and aggregates during concrete mixing, and conducted an analytical comparison between slump experiments and simulation to verify the feasibility of DEM method to study concrete materials.However, no systematic industrial application of the DEM method in and pipeline conveying has been seen.
Another analysis method is to use CFD-DEM (particle size small than grid).For example, Zhan et al. 21analyzed the trajectory of the particle motion in the pumping pipe and its exploration of the easy clogging parts of the pipe, and the results showed that the particles at the bends are prone to aggregation, and this result is consistent with the results analyzed by Jiang et al. 22 Considering the large grid size (large particle size), it is difficult to capture the action and characteristics of the boundary layer, which affects the final pressure calculation results.
In summary, the former research method of DEM is more for concrete materials and cannot characterize the fluid properties of concrete, so it is rarely involved in pumping aspects.The classical CFD method, which esteemed the concrete components as single-phase fluid, is far from the actual complex concrete composition.It is impossible to characterize the movement of large-grained aggregates in concrete.In addition, the interaction between aggregate and slurry cannot be characterized, and the causes of aggregate precipitation, settlement and blockage in the pumping process cannot be analyzed.It is necessary to couple the movement of slurry and aggregate.Therefore, the SPH-DEM coupling analysis method is selected here, which can consider the particle motion and the large rigid body motion of the concrete cylinder.
The smoothed particle hydrodynamics coupled with discrete elements method(SPH-DEM) is a sub-particle scale (SPS) analysis method, which can simulate large-scale particle motion (particle size larger than the grid) 23,24 compared to the classical computational fluid dynamics (CFD) method.The SPH-DEM method is used to establish a concrete model for the pumping rheological behavior of concrete, and the pumping experiments are used to test the applicability of the method and provide evaluation, prediction and improvement of analysis methods for pumping experimental results.

SPH method
Fresh concrete is a multi-component complex material.The scale can be divided into two phases, namely, largescale coarse aggregate solid phase and small-scale fluid phase.Smoothed particle hydrodynamics (SPH) [25][26][27] is used to solve the fluid phase, which is a computational fluid dynamics method under the Lagrangian framework.In SPH, 'particles' represent the fluid space domain.The solution of SPH is solved by the integral and differential of the flow field variables (such as velocity/pressure, etc.) and the kernel function.The Navier-Stokes (N-S) equations of computational fluid dynamics, which continuity and momentum equations are written as Eqs.(1) and (2)   where ρ is fluid density, p is fluid pressure, U = U_x, U_y, U_z is velocity vector, g = g_x, g_y, g_z is gravi- tational acceleration, Ŵ = µ∇ 2 U is viscous dissipation term.In SPH method, Ŵ can be characterized in two ways, artificial viscosity 28 and SPS turbulent viscosity 29 .
The solution of the Eqs.(1) and (2) consists of two main parts, which are the approximation of the kernel function and the approximation of the particle position, respectively.More detailed analysis of the kernel function can be found in the literature 30,31 .The integral of the function is transformed into the integral of the kernel function written as Where is the computational domain, W is the kernel function, which decreases monotonically with dis- tance, and h is the smooth kernel length, which determines the range of action of the kernel function, which corresponds to the discrete space as (1) www.nature.com/scientificreports/ in Eq. ( 4), m j is particle mass and ρ j is density.Similarly, the differentiation of the function can be written as where ∇ i is the difference value calculated with respect to particle i position.The SPH kernel function used in this study is the Wendland Quintic fifth-order interpolation kernel function where, α D is a numerical value with respect to dimensionality and is taken as 7/4πh 2 for the 2-D case and 21/16 π/h 3 for the 3-D case. the image of the kernel function and its corresponding differential counterpart are shown in Fig. 1.Therefore, the N-S equation in SPH scheme 26 is written as where, p a is the pressure at particle a and ∇ a denotes the gradient with respect to the coordinates of particle a. γ = 7, b = c 2 0 ρ 0 /γ , ρ 0 is the reference density and which is the speed of sound at the reference density.
Considering that when the fluid is in contact with a solid or a boundary, the integration region of the kernel function is not the whole region that the fluid can contain at this time, the range of action of the kernel function will be truncated by the boundary, which will lead to the loss of accuracy of the kernel function of the SPH 31 .To solve this problem the Dynamic boundary condition(DBC) is chosen for the work of this paper, while DBC is used in the study of Newtonian and non-Newtonian fluids 14,32,33 .
DualSPHysics (open-source CFD code based on SPH method) v5.0.51 [34][35][36] is used to solve computational fluid dynamics.The Project Chrono engine (open-source DEM code) is coupled for collisions and interactions between solid large particles and bottom area, rigid bodies, and other particles.

DEM
There are a large number of coarse aggregate particles inside the fresh concrete material.The motion process of large particles is generally solved by DEM method, and then the velocity and displacement are solved by Newton's law of motion.The solution process of the DEM method is shown in Fig. 2. DEM discretizes the large particles of coarse aggregate in concrete into many small spheres, and obtains the force of the coarse aggregate movement process by solving the collision and deformation between the spheres.
In Fig. 2, the force on the particles can be decomposed into tangential force and normal force.The two forces are composed of the elastic collision force and the damping dissipation force of the material, and the normal force can be obtained from the Coulomb friction law.
(5) where, e ij is the average elastic collision coefficient of the two materials.
where R ij is the particle radius, E ij is the Young 's modulus of the material, and v p1 is the Poisson 's ratio of the material.For the simplification of the model, the tangential force of large solid particles is taken as a certain proportion of the normal force.

Coupled SPH-DEM method
It is also possible to derive the movement of an object by considering its interaction with fluid particles and using these forces to drive its motion.This can be achieved by summing the force contributions for an entire body.By assuming that the body is rigid, the net force on each boundary particle is computed according to the sum of the contributions of all surrounding fluid particles according to the designated kernel function and smoothing length.Each boundary particle k therefore experiences a force per unit mass given by Eq. ( 15) where, f ka is the force per unit mass exerted by the fluid particle a on the boundary particle k , which is given by For the motion of the moving body, the basic equations of rigid body dynamics can then be used where M is the mass of the object, I the moment of inertia, V the velocity, the rotational velocity and R 0 the centre of mass.Equations Eqs. ( 17) and ( 18) are integrated in time in order to predict the values of V and for the beginning of the next time step.Each boundary particle within the body then has a velocity given by ( 11) www.nature.com/scientificreports/

Rheological model of concrete
Bingham and H-B models were used to study the rheological properties of fresh concrete.Both of them consider the initial flow state of concrete, and the shear stress τ should be greater than the yield stress τ y of concrete.The mathematical expression of Bingham model 37 and H-B model is where,τ B is the shear stress under the Bingham model,τ y is the yield stress, µ p is the plastic viscosity, γ is the shear rate, τ H−B is the shear stress under the H-B model, K( γ ) is the consistency coefficient, and n is the fitted power law index.In Eqs. ( 14) and ( 15), the part of shear stress less than yield stress is defined by viscosity as the viscosity of the motion process is When the shear rate is small, the viscosity tends to infinity at this time, which will bring convergence problems to the simulation process.Therefore, it is necessary to improve the Eq. ( 15) reasonably to make the model converge.The improved method used in this paper is to refer to the idea proposed by Papanastious et al. 38 , and make a reasonable correction to convert the H-B model from a discontinuous process to a continuous HBP process, such as The viscosity of the motion process at this time is where, Therefore, the rheological model under the HBP model shows continuity, and as the management parameter m increases, it will become more and more close to the ideal Bingham model, and its mapping relationship and image are shown in Fig. 3.
The Bi-viscosity model is based on the HBP model, which enlarges the viscosity coefficient in the small shear rate region, making the Bi-viscosity model closer to the ideal Bingham model.The image is shown in Fig. 4. It (20)

Stress analysis of concrete rheometer
The force diagram of the concrete rheometer during the sampling process is shown in Fig. 5.The blade size of the rheometer is fixed, the height of the cylinder is fixed, and the height of the support frame is fixed.It is assumed that the geometric size is AB = D , IL = D t , FI = z 2 , EF = h , ME = z 1 .The areas driven by the rheometer are divided into six areas: ABCD, EFAC-BGHD, ABON, CDKJ, FCJI-DHLK, and MNAE-OPBG.The torques of the above six regions are calculated respectively.The torque of the rotary process of the rheometer is where, T 1 represents the reaction torque of concrete to the blade in the ABCD region, T 2 is the reaction torque of the EFAC-BGHD annular region to the rheometer, T 3 is the CDKJ torque in the top region, T 4 is the ABON torque in the bottom region, T 5 is the torque received by the bottom annular region of the FCJI-DHLK region, and T 6 is the torque received by the top MNAE-OPBG region.
Since the tangential velocity of the rotation axis is zero, the tangential velocities of the blades AC and BD are V = ωr , where ω(ω = 2πn ) represents the angular velocity of the blade rotation.The corresponding shear rate is defined as www.nature.com/scientificreports/Torque of regional ABCD The torque of the blade AB and CD is Previous studies [39][40][41] have shown that no shear behavior occurs in the ABCD region when the blade is rotating.When the blade rotates, the surface of the blade exerts a tangential force along the surface direction on the concrete in the ABCD region.The centrifugal force generated by the rotating fluid in the ABDC region is completely balanced by the fluid outside the ABDC region, and there is no radial displacement.After the rotational speed of the rheometer shaft is stable, the rotational speed of the fluid in the ABDC region is consistent with the rotational speed of the blade, and the fluid and the blade have no relative movement trend in the radial and tangential directions.Therefore, the fluid in the ABDC region does not apply shear stress to the blade, that is, the required torque T 1 = 0.

Ring region EFAC-BGHD torque T 2
This region can be simplified as an annular flow, and the speed of the blade is AC/BD.The speed of the blade is the boundary speed V = ωr , and the effective width of the region is e = (D t − D)/2 .At this time, the shear rate γ of the region is V (e) , so the corresponding torque T 2 is Torque T 3 of region CDKJ This area can be simplified into two rotating plates, and the upper part of the plate rotates unilaterally, for the occurrence of concrete shear.The velocity V along the radius direction on the plate is ωr , so it's shear rate γ is The torque integral T 3 acting on the micro-element in this region is The force analysis of the torque in this region is similar to that of T 3 , and the integral of this region is Torque T 5 of region FCJI-DHLK

Considering that the velocity of any height z in this region is
The corresponding shear rate is defined by the shear rate as The torque contribute in this region is Torque T 6 of region MNAE-OBGP The analysis process of this region is similar to T 5 , so the value is www.nature.com/scientificreports/0) and (31) with T 2 = T 3 and T 5 = T 6 , and T 1 = 0 .Therefore, the total torque contributed by the test area to the rheometer is When n = 1, the fitting method of Eq. ( 33) corresponds to the Bingham model, which simplifies the relation- ship between torque and speed, that is, The calculation result of Eq. ( 34) is the same as the Bingham fitting formula deduced by Laskar et al 42 , which indirectly proves the correctness of the derivation form.In the actual measurement process of the rheometer, it is not a linear relationship similar to the Bingham form, but a nonlinear Bi-viscosity form Eq. ( 33) is written as Substituting the geometric parameters of the rheometer blade h = 0.127 m, z = 0.114 m, D = 0.127 m, D t = 0.38 m to Eq. ( 36), the corresponding relationship between shear stress and torque can be obtained : The shear strain corresponding to rotation is Similarly, the shear rate of the Bingham model is

Experimental material
The materials used in the four groups of concrete experiments for this test are ordinary Portland cement PO42.5, with a density of 3050 kg/m 3 .The type of fly ash is F class two, and the density is 2240 kg/m 3 .The type of ore powder is S95, and the density is 2880 kg/m 3 .The type of coarse aggregate is natural pebble, and its continuous gradation is 4.75-20 mm.The average aggregate density measured by multiple drainage method is 2600 kg/m 3 .The fine aggregate is machine-made sand and stone, and the humidity is 3.5%.The measurement method is shown in Fig. 6.The additive is a high-efficiency retarding water reducer, and the parameter is 9% In order to comprehensively study the influence of the rheological properties of concrete on the flow process, different proportions of concrete were prepared in this paper to prevent contingencies during the experiment.The specific match is shown in Table 1.
Considering that the material used in the experimental process is less, it is necessary to accurately control the quality of the material, especially the control of the water content has a great influence on the experimental results.In addition, in order to simulate the movement of large aggregate in concrete, the size and gradation distribution of aggregate in concrete are screened as shown in Fig. 7.

Slump test
The slump test is to detect the pumpability of concrete.Generally, the fluidity of concrete with slump greater than 180 mm is higher.The specific structure of the slump instrument is a cone with a top diameter of 100 mm, a bottom diameter of 200 mm, and a height of 300 mm, as shown in Fig. 8a and b.
(37)     www.nature.com/scientificreports/ In order to ensure the relative stability of the data obtained during the collapse process, the above four experimental groups were stirred and tested as shown in Fig. 9, and the average value of the slump s and the two-direction expansion s f was taken.The average value of the two extensions is s af , and the experimental results are shown in Table 2.

Concrete rheological test
The torque of the rheometer blade during rotation is calculated to convert the rheological properties of the concrete.ICAR Rheometer is used, as shown in Fig. 10.During the test process, the rheometer is standing still, and the rheometer is sunk into the fresh concrete by gravity.When the rheometer is started, the initial speed of the blade will change to 0.55rev/s.After a period of operation, the speed will decrease by 0.075rev/s.After multiple decelerations, the minimum measurement speed is 0.15rev/s, and the torque value T of the rheometer at multiple test sampling points is recorded.
The rheometer in Fig. 10 is divided into several parts, which are laptop, rheometer, blade, support frame and container cylinder.The rheometer is assembled when used, and the blade is inserted into the fresh concrete as shown in Fig. 11a.The torque data of the rotation process is read and recorded, and the recording process is shown in Fig. 11b.
Considering that during the operation of the rheometer, if the blade stirs to a larger aggregate, a certain fluctuation will occur, which will affect the final calculation results and thus affect the analysis results.Therefore, it is necessary to process the data of the rheometer.The original data of the rheometer is shown in Fig. 12.   Obviously, the data in Fig. 12 fluctuate greatly when the rotational speed of the rheometer changes or the larger aggregate is touched.The processing method of the original data is to select the last 3 s after the rotational speed of the rheometer is stable as the torque at the current rotational speed.In addition, it is also necessary to filter out the data points with large fluctuations in the rotation process, and then take the average torque in the speed range.The torque distribution after filtering is shown in Fig. 13.
The rheological test data of the concrete in the experimental group are shown in Table 3. Substituting the data in the table into Eqs.(36), (38) and (39), the Bingham and Bi-viscosity model data corresponding to the experimental group can be obtained, as shown in Table 4.
The data in Table 4 are substituted into the Eqs.( 22), ( 23), ( 24), ( 26), (27), and the Bi-viscosity model of the experimental group is fitted as Table 5 and Fig. 14    www.nature.com/scientificreports/Similarly, the Bingham model fitting can be obtained as shown in Table 6, and the corresponding Matlab data fitting is shown in Fig. 15.

Simulation analysis of rheological behavior of concrete
The solver used for this work is the explicit integration solver Symplectic.The fluid-solid interaction part uses Inlet/outlet boundaries with an initial time step of 1e−6.The computation time step is automatically adjusted based on the CFL number.The GPU of the computer is NVIDIA RTX A60000, and the CPU has 96 cores.

Simulation analysis model of concrete
The simulation model of concrete slump is composed of two phases, which are slurry, generated by fluid particles of Bi-viscosity model, and large aggregate DEM solid.The aggregate file is scanned into stl format and converted into SPH particles to form rigid large particle Vtk format, as shown in Fig. 18.The numbers in the figure is the size of the aggregate (x-y-z) in three directions (mm).
The concrete slump instrument is filled with slurry and coarse aggregate, and the proportion of slurry and aggregate is calculated by the results given in Fig. 7. Since the position of the aggregate in slump test is random, the aggregate initial position in the simulation model also needs to be randomly generated.Here, the Matlab  Monte Carlo random method generates aggregate positions of various sizes in the slump instrument, as shown in Fig. 19.Simulation parameters shown in Table 7.

Simulation analysis of concrete slump
After the concrete slump model is established, the rheological parameters and aggregate friction coefficient of the experimental group are substituted for simulation analysis.The simulation and experimental results are compared.The first set of simulation analysis is shown in Fig. 20.In order to facilitate the distinction between coarse aggregate and slurry, slurry in Fig. 20 is represented by gray, and the other is shown by velocity legend.When the coarse aggregate stops moving, the slump height will hardly change at this time.Slurry velocity could of slurry in 5 s is shown in Fig. 21.
It can be seen from Fig. 22a that at 1 s, the movement speed of the lifting cylinder and the concrete will gradually decrease, from the initial high-speed movement, a sharp decline to a stable process.The simulation slump is 63 mm at 4 s, and only 1 mm at 5 s.Compared with the difference of 13 mm between 3 and 4 s, the change is very small, so the simulation time is 5 s.The simulation data and experimental data results are shown in Fig. 21b.It can be seen from the diagram that compared with the experimental group E 1 , E 2 , it can be seen that the Bi- viscosity model can achieve a corresponding high accuracy range similar to the Bingham model.Compared with the experimental values, the error ratio is less than 10%.Combined with the slump simulation results of Fig. 22b, it can be seen that the overall error ratio of the Bi-viscosity model does not exceed 10%.
The results of the third E 3 and fourth E 4 experimental groups are shown in Fig. 23a and b.Compared with the experimental groups E 3 and E 4 , in the case of low-speed motion, due to the characteristics of the Bingham model itself, a large yield stress is required to make the concrete easy to flow.Therefore, in this case, the Bingham model fitted by the HBP model has poor numerical fitting accuracy for the experimental group, and the relative error can reach 20%, while the Bi-viscosity model still has high fitting accuracy for the rheological experiment, and its error percentage can be less than 10%.The influence of rheological parameters on the slump results

Effect of yield stress on slump
There are three main parameters in the Bi-Viscosity model, which will affect the simulation results.Based on the Bi-viscosity rheological parameters of the experimental group E 1 the yield stress in the model was changed to simulate.Considering the Bingham model used in the empirical formula s = 255 − 176τ/ρ 44 , n takes 1 here for comparative analysis.Considering that the yield stress is related to the anti-deformation ability of the fluid, the greater the yield stress, the stronger the anti-deformation ability, and the slower the flow.The simulation analysis results are shown in Fig. 24.It can be seen that as the yield stress increases.
Comparing the simulation results with the empirical formula, it can be obtained that the simulated slumps with yield stress of 400, 200 and 100 are 230 mm, 231 mm and 237 mm, respectively, while the empirical formula results are 246 mm, 237 mm and 220 mm.Therefore, the error percentages of simulation analysis and empirical formula are 6.6%, 2.8% and 7.5%, respectively.The overall error percentage is less than 10%.

The influence of coarse aggregate shape on simulation results
The comparison between spherical coarse aggregate and non-spherical coarse aggregate is mainly to consider the interaction between aggregate shape and bottom and aggregate.Because the contact and action range of spherical coarse aggregate is less than that of non-spherical coarse aggregate.Therefore, the expansion range of spherical coarse aggregate should be larger than that of non-spherical aggregate.The simulation results are shown in Fig. 25.The results of Fig. 25 are consistent with the trend of theoretical analysis results.

The influence of power law coefficient on the simulation analysis of slump
The power law coefficient in the Bi-viscosity model is changed to explore the influence of its change trend on the concrete slump experiment, as shown in Fig. 26.It can be seen that when the initial acceleration changes greatly, the expansion range of the shear thinning fluid is obviously larger than that of the shear thickening fluid.With the extension of the simulation time, when the overall acceleration changes little, the corresponding shear stress is similar to the yield stress.Therefore, for the final result, it tends to the expansion of n = 1, and the overall impact is not obvious.

The influence of power law coefficient on the simulation analysis of slump
Changing the fluid density for slump simulation analysis, the most intuitive change brought by density is the difference in gravity, which leads to the difference in shear stress corresponding to the motion process.The greater the density, the greater the shear stress, and vice versa.Therefore, the slump and expansion of the fluid with higher density will increase, and the slump of the fluid with lower density is relatively small.The simulation results are shown in Fig. 27.

Comparison of simulation and experimental results of concrete pumping process
The industrial application of fresh concrete is mainly concrete pumping.In order to explore the applicability of the model in industrial application, a simple pumping platform is established, as shown in Fig. 28.The pumping pressure of the experimental group is compared with the pressure change in the simulation analysis of the pumping push process.If the matching degree is high, it can provide some solutions and methods for the coarse aggregate motion trajectory, pipe plugging reason, pressure prediction and analysis in the concrete pumping process.
by using Matlab data fitting tool.

Figure 14 .
Figure 14.Fitting curves of test data based on Bi-viscosity model.

Figure 15 . 6 Figure 16 .
Figure 15.Fitting curves of test data based on Bingham model.

Figure 17 .
Figure 17.Drag coefficient C d as a function of the Reynolds number Re.

Figure 20 .
Figure 20.Simulation result of first group E 1 .

Figure 21 .
Figure 21.Slurry velocity in Y and Z directions.

Figure 22 .
Figure 22.The velocity and spread change with time.

Figure 24 .
Figure 24.Effect of yield stress on slump test.

Figure 25 .
Figure 25.Comparison of spherical and non-spherical aggregate.

Figure 26 .
Figure 26.Influence of power-law coefficient on slump analysis.

Figure 27 .
Figure 27.Influence of fluid density on slump.

Figure 32 .
Figure 32.Comparsion of the simulation and experiment result.

Table 3 .
Relationship between torque and rotational speed in experimental group.

Table 4 .
Relationship between shear stress and shear rate.

Table 5 .
Rheological parameters of fresh concrete based on Bi-viscosity model.
y (Pa) K(Pa•s n ) n

Table 6 .
Rheological parameters of fresh concrete based on Bingham model.